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ABSTRACT 

Close interactions and mass transfer in binary stars can lead to the formation of many different 
exotic stellar populations, but detailed modeling of mass transfer is a computationally challenging 
problem. Here, we present an alternate Smoothed Particle Hydrodynamics approach to the modeling 
of mass transfer in binary systems that allows a better resolution of the flow of matter between 
main-sequence stars. Our approach consists of modeling only the outermost layers of the stars using 
appropriate boundary conditions and ghost particles. We arbitrarily set the radius of the boundary 
and find that our boundary treatment behaves physically and conserves energy well. In particular, 
when used with our binary relaxation procedure, our treatment of boundary conditions is also shown 
to evolve circular binaries properly for many orbits. The results of our first simulation of mass transfer 
are also discussed and used to assess the strengths and limitations of our method. We conclude that 
it is well suited for the modeling of interacting binary stars. The method presented here represents 
a convenient alternative to previous hydrodynamical techniques aimed at modeling mass transfer in 
binary systems since it can be used to model both the donor and the accretor while maintaining the 
density profiles taken from realistic stellar models. 

Subject headings: binaries: close — stars: evolution — hydrodynamics — methods: numerical 



1. INTRODUCTION 

Gravitation rules. It is what forms dark matter halos 
and giant molecular clouds. It is also what compresses 
these clouds of gas to such an extent that stars form 
out of them. And these newly-formed stars are often 
bound to each other by the means of this force, forming 
binary or multiple stellar systems. Although the details 
of bi nary star form ation are still not fully understood 
(e.g. iTohlind 120 02) . it is now acknowledged that stellar 
multiplicity is more the rule than the exception. Obser- 
vations suggest that over 50% of the stars in our Galaxy 
are part of double, triple, quadruple, or even sextuple 
systems (lAbt fc Levvi 119761: iDuauennov fc Mayor! Il991t 
"IT 



iFischer fc Marcvlll992t lHalbwachs et all 12003! ) ^Because 
stars grow in size considerably as they evolve, it is es- 
timated that those binaries with a period of less than 
~ 1 4 days will inevi tably interact at some point of their 
life (|Paczvriskilll971l ). When such close interactions oc- 
cur, material is transferred from one star to the other and 
the course of evolution of the stars is irreversibly altered. 
Likewise, such close interactions can also be triggered 
by stellar encounters in dense stellar environ ments (e.g. 
captures and exchanges: lPoolev""fc Hut 2006h. 

The first realization of th e importance of binary in- 
teractions may have been by I Crawford! ([1955), who sug- 
gested an interesting solution to the paradox of the Al- 
gol system, in which the more evolved star is also the 
least massive. ICrawforcH suggested that the close prox- 
imity of the two components must have led to significant 
mass transfer from the initially more massive star to the 
least massive one until the mass ratio was reversed. This 
discovery opened the way to many more types of stars, 
such as cataclysmic variables and X-ray binaries, helium 
white dwarfs, and blue stragglers, whose very existence 
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could now be understood in terms of close binary evo- 
lution. Even for stars that are not transferring mass, a 
close companion can have all sorts of effects on their ob- 
servable properties, such as increased chromosphcric and 
magnetic act ivity (e.g. RS Canum Venaticorum stars; 
lRodonolfl99l . 

Since the work of ICrawfordl ()1955l ) , much effort has 
been put into bet ter unde r stand i ng binary evoluti o n and 
its by-products (iMortonl Il960t iPaczvnskil 119651 119711: 
iPaczvnski fc SienkiewiczHl972tlIbenlll99lD . However, the 
usual Roche lobe formalism used to study binary stars 
applies only to the ideal case of circular and synchronized 
orbits. The addition of simple physics such as radiation 
pressure, for example, is enough to sig nificantly mod- 
ify th e equipotentials of binary systems (jDermine et al.1 
l2009f ) and estimates of fundamental parameters such 
as the Roche lobe radius Rl consequently become un- 
certain. Since the evolution of close binaries depends, 
among other things, on the rate at which mass transfer 
proceeds and Rl, analytical prescriptions for the mass 
transfer rate also become uncertain. The determination 
and characterization of the mass transfer rate in binary 
stars therefore represents a key issue that needs to be ad- 
dressed, especially given that surveys of binary stars have 
show n that a non-negligible fraction (~ 20% in the sam- 
ple oflPetrova fc Orlovlll99"9t see also lRaguzova fc Popov! 
2005]) of semi-detached or contact systems have eccentric 
orbits. 

Re cent analytical work by ISepinskv et al.l (|2007atf bl. 
l2009f ). who investigated the secular evolution of eccen- 
tric binaries under episodic mass transfer, have shown 
that, indeed, eccentric binaries can evolve quite differ- 
ently from circular ones. However, because mass transfer 
can occur on short dynamical timescales, it is necessary 
to use other techniques to characterize it fully. In partic- 
ular, hydrodynamics has shown to be useful for studying 
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transient phenomena and episod es of stable mass trans- 
fer. Simple ballistic models (e.g. IWarner k. PetersllT972l : 
iLubow fc Sh^[T^lFlannervl[l975ir and two-dimensional 
hydr odynamical simulat i ons of semi-detache d binaries 
(e.g. ISawada et aTlll986l iBlondin et aTlll995l ) have gen- 
erally been used in the past to study the general char- 
acteristics of the flow between two stars and the proper- 
ties of accretion disks. Later three-dimensional models 
with higher resolution allowed for more realistic stud- 
ies o f coalescing binaries (e.g. iRasio fc Shapiro] 11994 
1 1995T) and accretion disks (e.g. iBisikalo et al.l 119981 ) in 
semi-detached binaries, all mainly focusing on the sec- 
ular and hydrodynamical stability of binaries and on 
the structure of the mass tra nsfer flow. More recently , 
iMotl. Tohline. fc Frankl d200l and lD'SouzaeTall $MM 
used grid-based hydrodynamics to simulate the coales- 
cence of n = 3 /2-poly tropes, representative of low- mass 
main-sequence stars (M < 0.5 M Q ). They also investi- 
gated the onset of dynamical and secular instabilities in 
close binaries and were able to get a good agreement with 
theoretical expectations. 

Characterization of the accretion process and the be- 
haviour of the accreted material onto the secondary star 
requires that the secondary be modeled realistically. In 
most simulations to date, the secondary has often been 
modeled using point masses or simple boundary condi- 
tions to appro ximate the surface accret ion. Moreover, as 
pointed out bv lSills &: Lombard! (|1997t ). the use of poly- 
tropes, instead of realistic models, may lead to signifi- 
cantly different internal structures for collision products, 
which may arguably be applicable to mass-transferring 
binaries. Therefore, more work remains to be done in 
order to better understand how mass transfer operates 
in binary systems. 

Here, we present our alternate hydrodynamical ap- 
proach to the modeling of mass transfer in binaries. We 
first discuss the usual assumptions made when study- 
ing binary stars 3D Our computational method, along 
with our innovative treatment of boundary conditions, 
are introduced and tested in $3] We discuss initial con- 
ditions for binary stars in SPH and the applicability of 
our technique to binary stars in fj4] Our first mass trans- 
fer simulation is analyzed in $5l The work presented 
here is aimed at better understanding the process of mass 
transfer and wi ll be applied to spec ific binary systems i n 
another paper (|Laioie fc Sillsll2010l; hereafter, |PaperL|). 

2. BRIEF THEORY OF BINARY SYSTEMS 

Our understanding of interacting binary stars is 
prompted by observations of systems whose existence 
must be a result of mass transfer. These observations 
and simple physical considerations have driven the theo- 
retical framework we use to model these systems. 

2.1. The Roche approximation 

Under the assumptions of point masses (Mi and M2), 
circular and synchronous orbits, the gravitational poten- 
tial in a rotating reference frame gives rise t o equipo- 
tentia ls with particular shapes, as described in lEggletonl 
(|2006f l. The equipotential that surrounds both stars and 
intersects at a point between the two stars is the Roche 
lobe and forms a surface within which a star's potential 
dominates over its companion's. lEggletonl (|1983l ) approx- 



imates the equivalent Roche lobe radius (Rl) to an ac- 
curacy of 1% over the whole range < q± < 00, by 
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where a is the semi-major axis. Analytical studies of 
binary stars usually rely on this definition of the Roche 
lobe radius, but as discussed in [JU such estimates for Rl 
may not always be reliable. 

2.2. Rate of mass transfer 

Once a star overfills its Roche lobe, mass is transferred 
to the companion's potential well. The rate at which 
material is transferred is, in general, a strong function of 
the degree of overflow of the donor, AR = Rl, where 
i?* is the radius of the star. In general, mass is assumed 
to be transferred through the Li point and, using simple 
physical assumptions such as an iso therma l and i nviscid 
flow and an ideal gas pressure law, iRitterl (|1988l ) shows 
that the mass transfer rate can be expressed as 



Mi = -M exp 



^ph _ r l 



(2) 



Here, R p h is the photosphcric radius, Hp is the pres- 
sure scale height of the donor star, and Mo is the mass 
transfer rate of a star exactly filling its Roche lobe. This 
mass transfer rate depends rather strongly on the de- 
gree of overflow AR and goes to z ero exp o nenti ally if 
the star is within its Roche lobe. Rittcr (1988) pro- 
vides estimates for Mq of about 10 -8 M yr" 1 and 
Hp/R ~ 10 -4 for low-mass main-sequence stars. This 
model of mass transfer has been successfully applied to 
cataclysmic variables where the photosphere of the donor 
is located about one to a few Hp inside the Roche lobe. 
For cases where the mass trans fer rates a re much larger, 
iPaczvriski fc Sienkiewici (Il972ft (see also lEggletonl 120061 
and Gokhale et al. n2007f) derives, in a similar way, the 



mass transfer rate for donors that can be approximated 
by polytropes of index n. In such cases, the mass transfer 
rate is 

• /i?^ — Rt \ n+3/2 
Mi = -M ( U *„ Ul ) (3) 



where Mq is a canonical mass transfer rate which depends 
on Mi, M2, and a. The dependence of the mass transfer 
rate on (AR/R) is again found, although somewhat dif- 
ferent than that of IRitterl ((l988(). This rate is also zero 
when AR < and is applicable when the degree of over- 
flow is much larger than the pressure scale height and 
mass transfer occurs on a dynamical timescale. 

In any case, once a star fills its Roche lobe, the mass 
transfer is driven by the response of the star's and Roche 
lobe radii upon mass loss. Stars with deep convective 
envelopes tend to expand upon mass transfer whereas 
radiative stars tend to shrink upon mass transfer. The 
ensuing response of the Roche lobe radius, which may 
expand or shrink, will therefore dictate the behaviour of 
the mass transfer rate. 

2.3. Orbital evolution due to mass transfer 

If the mass of the stars changes, then both the period 
and separation ought to readjust. This behaviour can be 
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shown by taking the time derivative of the total angular 
momentum of a system of two point mass orbiting each 
other with an eccentricity e: 



Ltot 



Mi M 2 1 M , la 



ee 



(4) 

Mi ' M 2 2 M ' 2 a 1 - e 2 " w 

Equation 2] shows that as the masses of the stars change 
and as mass and angular momentum are being lost from 
the system, both the orbital separation and the eccen- 
tricity change. The exact behaviour of these quantities 
depends of course on the degree of conservation of both 
total mass and angular momentum. For the usual as- 
sumptions of circular orbits and conservative mass trans- 
fer, we can further impose that e = 0, Mi — —M\ and 
•hot = 0, therefore reducing Equation [4] to 



a Mi V Mi 



(5) 



Assuming Mi to be the donor and more massive (i.e. 
Mi < and Mi > M<x) we find that the separation a de- 
creases until the mass ratio is reversed, at which point the 
separation starts increasing again. For main-sequence bi- 
naries, where the most massive star is expected to overfill 
its Roche lobe first, the separation is therefore expected 
to decrease upon mass transfer. 

2.4. Limitations 

The theoretical framework derived in this section has 
generally been applied to the study of close binaries. 
However, strictly speaking, it is not valid in most in- 
stances. Close binaries are not all circular and synchro- 
nized, and the Roche lobe formalism therefore does not 
apply. This, in turn, makes estimates of mass trans- 
fer rates rather uncertain. Moreover, conservative mass 
transfer is more an ideal study case than a realistic one 
and the secular evolution of binary system becomes a 
complex problem. To circumvent these difficulties, ap- 
proximations to the mass transfer and accretion rates as 
well as to the degree of mass loss have to be made. How- 
ever, to better constrain these approximation or avoid 
over-simplifications, one can use hydrodynamics, which 
is well suited for modeling and characterizing episodes of 
mass transfer. We therefore discuss our hydrodynamics 
technique and show how it can be used to better con- 
strain mass transfer rates in binary systems. 

3. COMPUTATIONAL METHOD 

3.1. Smoothed Particle Hydrodynamics 

Smoot hed Particle Hydr odynamics (SPH) was intro- 
duced bv lLucvl (| 1977ft and iGingold fc Monaghanl (|1977l ) 
in the context of stellar astrophysics. Its relatively sim- 
ple construction and versatility have allowed for the 
modeling of ma ny different physica l problems such a s 
star formation dPrice fc Bate! l2009t iBate et al.1 I1995D . 



accretion disks 



stellar col- 



.Maver et al.l .. 
lisions (lLombardi et al.l 119951 iSills et al.1 fl997l 12001 . 

galaxy formation and cosmological simulations (e.g. 



Mashchenko. Couchman. fc Wadslcv 2006; IStinson et all 
2009trG overna to et al.l [2009). Our code derives from that 
of iBatd (11995f) . which is bas e d on the earlier version of 
iBenzl (|1990¥ and lBenz et al.1 (|1990H . Here, we only em- 
phasize on the main constituents of our code. The reader 



is referred to these early works for complementary de- 
tails. 

SPH relies on the basic assumption that the value of 
any smooth function at any point in space can be ob- 
tained by averaging over the known values of the func- 
tion around this point. This averaging is done using a 
so-called 'smoothing kernel' to determine the contribu- 
tion from neighbouring particles. The smoo thing ker- 
nel can take many forms (see e.g. iPricd 120051 ): here we 
use the com pact and spherically symmetri c kernel first 
suggested bv lMonaghan fc Lattanziol ([1985D . To prevent 
the intcrpenctration of particles in shocks and allow for 
the dissipation of kinetic energy into heat, we include 
an artificial viscosity term in the momentum and energy 
equations. The artificial viscosit y can also take various 
forms (e.g.lLqmbardi et al.l fl999h: we use the form given 
bv lMonaghan! (|1989t ) witha = 1 and (3 = 2. We allow for 
the smoothing length to change both in time and space, 
and we use individual timesteps for the evolution of all 
the required quantities. In this work, we assume an equa- 
tion of state for ideal gases of the form P = (7 — l)pu, 
where 7 = 5/3 is the ratio of the heat capacities. Finally, 
we use the parallelized version of our code (OpcnMP), 
which scales linearly up to ~ 24 CPUs for simulations of 
~ 10 6 particles. 

3.2. Boundary conditions 

As discussed bv lDeupree fc Karakasl (|2005f ). the inner 
parts of stars in close binaries generally remain unaf- 
fected by the presence of a companion, and only the 
structure of the outermost layers is modified by close 
tidal interactions. This result prompted us to model 
only the outer parts of the stars with appropriate bound- 
ary conditions. Such an approach effectively reduces the 
total number of SPH particles in our simulations with- 
out decreasing the spatial resolution. Conversely, for the 
same amount of CPU time, modeling only the outermost 
layers of stars allows for the use of more particles, there- 
fore enhancing the spatial and mass resolutions. More- 
over, CPU time is spent solely on particles actually tak- 
ing part in the mass transfer or being affected by the 
companion's tidal field. 

SPH codes calculate hydrodynamical quantities by av- 
eraging over a sufficiently large number of neighbours. 
For particles located close to an edge or a boundary, two 
things happen. First, since there are no particles on one 
side of the boundary, a pressure gradient exists and the 
particles tend to be pushed further out of the domain 
of interest. Second, if the number of neighbours for each 
particle is kept fixed by requirements, as it is in our code, 
then the smoothing length is changed until enough neigh- 
bours are enclosed by the particle's smoothed volume. 
This lack of neighbours therefore effectively decreases the 
spatial resolution at the boundary and underestimates 
the particle's density. In such circumstances, the imple- 
mentation of boundary conditions is required. 

3.2.1. Ghost particles 

Boundary conditions have often been implemented us- 
ing the so-called ghost particles , first intr oduced by 
Takeda. Mivama. fc Sekival (|1994[ ) (see also iMonaghanl 
1994f ). Ghost particles, like SPH particles, contribute 



to the density of SPH particles and provide a pressure 
gradient which prevents the latter from approaching or 
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FIXED GHOSTS 



point mass (since the point mass and the ghosts move 
together), which we write as 



Fig. 1. — Illustration showing how boundary conditions are im- 
plemented using fixed ghost particles. SPH and ghost particles are 
represented by solid and open dots respectively and the boundary 
is represented by the thick solid line. The density calculation of 
SPH particles includes the contributions from all particles located 
within two smoothing lengths of them. 



penetrating the boundary. Ghost particles can be cre- 
ated dynamically every time an SPH particle gets within 
two smoothing lengths of the boundary. When this oc- 
curs, the position of each ghost is mirrored across the 
boundary from that of its parent SPH particle (along 
with its mass and density). Therefore, the need for 
ghosts (and a boundary) occurs only when a particle 
comes within reach of the boundary. Here, however, 
w e use a slightly different ap proa ch, based on the work 
of iMorris. Fox, fc Zhul () 19971 ) and ICummins fc Rudmanl 
(|1999f l. The approach of these authors differs from the 
mirrored-ghost technique in that the ghosts are created 
once, at the beginning of the simulation, and their rela- 
tive position remains fixed in time during the simulation. 
We further improve upon this technique in order to model 
the outer parts of self-gravitating objects. Starting from 
our relaxed configurations, we identify any particles as 
ghosts if they are located within three smoothing lengths 
inside of the boundary, which, at this point, is arbitrar- 
ily determined. Particles located above the boundary 
are tagged as SPH particles, whereas the remaining ones 
are erased and replaced by a central point mass whose 
total mass accounts for both the particles removed and 
those tagged as ghosts. Point masses interact with each 
other and SPH particles via the gravitational force only. 
Point masses are also used when modeling massive or gi- 
ant branch stars whose steep density profile in the core 
is hard to resolve with SPH particles. Figure [T] illus- 
trates how our boundary conditions are treated in our 
code. Here, we use three smoothing lengths of ghosts 
as a first safety check in order to prevent SPH particles 
from penetrating the boundary. We also enforce that no 
particle goes further than one smoothing length inside 
the boundary by repositioning any such particle above 
the boundary. 

We further ensure conservation of momentum by im- 
parting an equal and opposite acceleration to the central 



a, 



E 



(6) 



where a^ is the hydrodynamical acceleration imparted 
to particle i from ghost j. This term is added to the 
usual gravitational acceleration of the central point mass. 
Ghosts are moved with the central point mass and are 
given a fixed angular velocity. Note that at this point, aj 
has already been calculated by our code so that this cal- 
culation requires no extra CPU time. Ghost particles are 
also included in the viscosity calculations to realistically 
mimic the interface. 

3.2.2. Applications to single stars 

We now show that our new boundary condition treat- 
ment is well suited for the modeling of stars in hydro- 
static equilibrium. We first relax a O.8-M0 star with 
rotation (w = 0.10; solar units). The relaxation of a star 
requires a fine balance between the hydrodynamical and 
gravitational forces, and therefore allows to assess the ac- 
curacy of our code. We model our stars using theoretical 
density profile s as given by the Yale R otational Evolution 
Code fYREC: [Guenther et al.lfl992T) . SPH particles are 
first spaced equally on a hexagonal close-packed lattice 
extending out to the radius of the star. The theoretical 
density profile is then matched by iteratively assigning a 
mass to each particle. Typically, particles at the centre of 
the stars are more massive than those located in the outer 
regions, by a factor depending on the steepness of th e 
density profile. As discussed bv lLombardi et all (|1995| ). 
this initial hexagonal configuration is stable against per- 
turbations and also tends to arise naturally during the 
relaxation of particles. Stars are relaxed in our code 
for a few dynamical times to allow for the configura- 
tion to redistribute some of its thermal energy and settle 
down. Once the star has reached equilibrium, we remove 
the particles in the central regions and implement our 
boundary conditions. We give the star, the ghosts and 
the central point mass translational and angular veloci- 
ties. The final relaxed configuration of our star, with the 
boundary set to ~ 75% of the star's radius, is shown in 
Figure [2] along with the density and pressure profiles in 
Figure [3] By using our initial configuration for the setup 
of ghosts, we ensure that the ghosts' position, internal 
energy, and mass are scaled to the right values and that 
the ensuing pressure gradient maintains the global hy- 
drostatic equilibrium. The total energy (which includes 
the gravitational, kinetic, and thermal energies) for the 
model of Figure [2] evolved under translation and rotation 
is conserved to better than 1% over the course of one full 
rotation. These results suggest that our treatment of 
boundaries is adequate for isolated stars in translation 
and rotation. 

4. APPLICATIONS TO BINARY STARS 

We now discuss the modeling of binary stars with our 
new boundary conditions. In particular, we develop a 
self-consistent technique for relaxing binary stars and 
show that our code, along with our new boundary con- 
ditions, can accurately follow and maintain two stars on 
relatively tight orbit for many tens of dynamical times. 
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Fig. 2. — Example of a star modeled with ghost particles (small 
red dots), SPH particles (black dots), which are found beyond the 
boundary (green line), and a central point mass (big red dot). 




Fig. 3. — Density and pressure profiles of the ghosts (red) and 
SPH particles (black) showing the gradients at the boundary, as 
expected for realistic models of stars. In all panels, only the par- 
ticles located within 2h of the equatorial plane are plotted. The 
profiles of the ghost particles exactly match that of the theoreti- 
cal profiles, showing that the overall profiles are very close the the 
theoretical ones. 



We emphasize that the location of the boundary is, at 
this point, arbitrary. We will discuss this issue in more 
details in $L3l 



4.1. Binary star relaxation 

As discussed in stars must be relaxed prior to 

being used in simulations. This is also true for binary sys- 
tems since tidal effects are not taken into account when 
calculating the theoretical density profiles of the individ- 
ual stars. Therefore, care must be taken when preparing 
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Fig. 4. — Logarithm of density in the orbital plane for a de- 
tached relaxed binary with stars of 0.80 Mq fully modeled (left) 
and modeled with boundary conditions (right). 



binary systems for SPH simulations. 

To account for the different hydrostatic equilibrium 
configuration of binary stars, we use the fact that the 
stars should be at rest in a reference frame that is cen- 
tered at the centre of mass of the binary and that rotates 
with the same angular velocity as the stars. A centrifugal 
term is added to the acceleration of the SPH particles, 
which we find by requiring that it ca ncels the net gravita- 
tional acceleration o f the stars (e.g. iRosswog et al.ll2004l : 
iGaburov et all2010l) . By assuming an initial orbital sep- 
aration, we find, at each timestep, the necessary angular 
velocity Q that cancels the stars' net gravitational accel- 
eration. Note that we do not account for the Coriolis 
force since the system is assumed to be at rest in the ro- 
tating frame. The net acceleration of the centre of mass 
of star i is calculated in the following way: 



hyd 



grav 



Mi 



(7) 



where Mi, a.h y d and a. grav are respectively the total mass 
and the hydrodynamical and gravitational accelerations 
of star i, and the summation is done over particles j 
that are bound to star i. We therefore get the following 
condition for the angular velocity: 



(8) 



3 3 



where rj is the distance of particle j to the axis of ro- 
tation of the binary. We find the angular velocity for 
both stars and take the average value and add it to the 
acceleration of each SPH particle. To ensure that the 
orbital separation is kept constant during the relaxation, 
we also reset the stars to their initial positions after each 
timestep by a simple translation. An example of such 
a detached relaxed binary is shown in the left panel of 
Figure [U 

4.2. Circular orbits 

After the initial relaxation, the stars are put in an in- 
ertial reference frame by using the value of the angular 
velocity obtained from the binary relaxation procedure. 
To assess the physicality and numerical integration ca- 
pabilities of our hydrodynamical code, we have evolved 
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Fig. 5. — Normalized energies and orbital separation as a function 
of time for a 0.6+0.6 Mq circular binary modeled with ~200,000 
particles initially relaxed using the method of £|4,ll with fixed orbital 
separation. 



different wide binaries on circular orbits for a small num- 
ber of orbits. This is important for simulations of mass 
transfer as any changes in orbital separation must be 
driven by the mass transfer itself and not the initial con- 
ditions and/or numerical errors inherent to our method. 

Figure [3] shows the normalized orbital separation and 
different energies as a function of time for a 0.60 + 0.60 
M binary with each star fully modeled (i.e. without a 
boundary). Each star contains ~ 105,000 SPH particles 
and the initial separation is 3.25 Hq. Our relaxation pro- 
cedure yields an angular velocity of il = 0.1875 which, 
because of the large initial separation, is identical to the 
Keplerian value. The different forms of energy are all 
well conserved, as well as the orbital separation, which 
remains constant to better than ~ 0.25% for over six 
orbits. The slight (anti-symmetric) variations observed 
in the orbital separation and the kinetic energy are an 
indication that the system is on a slightly eccentric or- 
bit. At this level, however, we estimate that our code 
can properly (and physically) evolve two stars orbiting 
around each other. 

Figure [6] shows the comparison between the normal- 
ized energies and orbital separation of a 0.80 + 0.80 M Q 
binary modeled with a different number of particles (the 
high resolution simulation is also shown in Figure EJ. 
The solid and dotted lines correspond, respectively, to 
the systems modeled without and with boundary con- 
ditions. In the simulations shown here, the boundary 
is set at 75% of the star's radius. Figure [5] shows that 
introducing the boundary conditions smoothes the os- 
cillations seen for the fully modeled binaries. This is 
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Fig. 6. — Normalized total energy and orbital separation for a 
0.80+0.80 Mq circular binary modeled with (dotted line) and with- 
out (solid line) boundary conditions. The upper two panels arc for 
a low-resolution simulation containing ~ 40, 000 (full star) parti- 
cles and the lower two panels are for a high-resolution simulation 
containing ~ 240,000 particles (full star). 



explained by the fact that the calculation of the gravita- 
tional force on the central point mass, and hence most 
of the star's mass, is done using a direct summation (in- 
stead of a binary tree) , therefore improving the accuracy 
and reducing the oscillations in the orbital separation. 
Figure [5] also shows that increasing the spatial resolution 
increases the quality of the evolution of the orbits. For 
comparison, the orbital separation of our high-resolution 
simulation varies by less than 0.5%, compared to ~ 2% 
for the low-resolution (Figure [S]). Moreover, in all cases, 
the total energy is conserved to ~ 0.25% over the whole 
duration of the simulations and the gravitational and 
thermal energies (not shown) remain constant to better 
than - 0.5%. 

Similarly, Figure [7] shows the evolution of a 0.80 + 0.48 
M.q binary system over four orbits and a relatively low 
number of particles. Note that only the primary is mod- 
eled with the use of our boundary conditions for rea- 
sons that are discussed in |J6J Again, the orbital sepa- 
ration remains fairly close to its initial value, to within 
2% at the end of four orbits. The total energy is con- 
served to better than 0.5% and only the kinetic en- 
ergy oscillates significantly. We note that, when com- 
paring our results of circular orbits with those of other 
authors, our binary relaxation procedure yields quan- 
tit atively comparabl e orbital behaviours. The results 
of iBenz et al.l (|1990t ) show oscillations of ~ 1 — 2% in 
the orbital separation over three orbits whereas the sim- 
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Fig. 7. — Comparison of the normalized energies and orbital 
separation as a function of time for a 0.80+0.48 Mq circular bi- 
nary modeled with (dotted line) and without (solid line) boundary 
conditions. 



ulations from iDan. Rosswog. fe Briiggcn (|2008| ) of two 
unequal-mass binaries showed a constant orbital sep- 
aratio n for many tens of orbits wi th a n accuracy of 
~ 1%. iMotl. Tohline. fc Franl d200l ) and lD'Souza et all 
using a specifically designed grid-based hydrody- 
namics code, both maintain equal- and unequal-mass bi- 
naries on circular orbits with an accuracy of ~ 0.25% 
over five orbits. 

Finally, we note that although our boundary condition 
treatment does require more calculations (e.g. the contri- 
bution from all the particles to the point mass' accelera- 
tion), the simulations of binary systems can be sped up 
by up ~ 50% when modeled with boundary conditions. 
This depends of course on the location of the boundary, 
which at this point is arbitrary but chosen wherever it 
makes the most sense for the problem at hand. Note, 
however, that the boundary should be placed at least a 
few smoothing lengths from the surface of the star. 

4.3. Eccentric Orbits 

The case of an eccentric orbit is interesting since tidal 
forces are time-dependent and can vary greatly depend- 
ing on the orbital phase. Here, we show that despite 
the large tidal forces at periastron, our boundary con- 
ditions are well suited for the modeling of such systems. 
The system we model consists of two main-sequence stars 
with masses of 1.40 and 1.50 M Q with an eccentricity of 
e = 0.15 and evolved for over four orbits (-P orD — 44 code 
units). The total number of particles nears ~ 500,000 
and the location of the boundary is at ~ 75% of the stars' 
radius, which, as shown in Figure [U is deep inside the 
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Fig. 8. — Radii enclosing different fractions of the total bound 
mass (in SPH particles) to the primary star as a function of time 
for the 1.40 + 1.50 M binary with e = 0.15. The dotted line 
represents the location of the boundary. 



star so that the effects of tidal force are negligible. In- 
deed, Figure |5] shows the radius of the primary enclosing 
different fractions of the total bound mass (in SPH par- 
ticles) as a function of time for our binary system. For 
example, the radii containing 60% to 90% of the total 
bound mass are shown to not change significantly dur- 
ing the whole duration of this simulation. In fact, only 
the outer radius of the star, containing over 95% of the 
bound mass, oscillates during each orbit. Therefore, in 
this case, the choice of the location of the boundary (dot- 
ted line) is well justified and Figure [8] shows that the use 
of our method for eccentric binaries is adequate. Replac- 
ing the core of a star with a central point mass and a 
boundary remains a valid approximation as long as the 
boundary is deep enough inside the envelope of the star. 



5. SIMULATION OF MASS TRANSFER 

We now present the results from the simulation pre- 
sented in i j4.3l In particular, we are interested in the 
mass transfer rates observed along the eccentric orbit. 
We also assess the physicality and limits of our approach 
later in ^ 

Figure O shows the logarithm of the density for parti- 
cles close to the orbital plane. The use of our boundary 
conditions can be seen as the centre of the stars is de- 
void of particles, except for the central point masses (not 
shown). Short episodes of mass transfer are observed 
close to periastron while mass transfer stops as the stars 
get further apart. The material being transferred hits 
the secondary and disturbs its outer envelope such that 
the latter loses some material. In the end, the secondary 
is surrounded by a relatively thick envelope. 

5.1. Determination of bound mass 

To estimate the mass transfer rates, we have to de- 
termine the component to which every SPH particle is 
bound. We use a total-energy (p e r unit mass) criterion, 
as presented by lLombardi et al.l (|2006l ) , and determine 
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Fig. 9. — Logarithm of the density in the orbital plane for the 
1.40 + 1.50 Mq binary with e = 0.15. The orbital period is about 
44 time units and the central point masses are not shown. 



whether a particle is bound to the primary, the sec- 
ondary, or the binary as a whole. In particular, given 
that most of the mass of the two stars is contained in 
the point masses, we use the latter as the main compo- 
nents to which particles are bound. For a particle to be 
bound to any of the components, we require its total en- 
ergy relative to the component under consideration to be 
negative. The total energy of any particle with respect 
to both the point masses and the binary's centre of mass 
is defined as 



E-, 



1 2 

2% 



G(Mj 



(9) 



where Vij and dij are the relative velocity and separa- 
tion, respectively, between particle i and component j. 
Moreover, we require the separation dij to be less than 
the current separation of the two centres of mass of the 
stars (in this case, the point masses). For particles that 
satisfy both of these criteria for both stellar components, 
we assign them to the stellar component for which the to- 
tal energy is most negative. If only the energy condition 
is satisfied for the stellar components, or the energy with 
respect to the binary is negative, the particle is assigned 
to the binary component. Finally, if the total relative en- 
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Fig. 10. — Change in bound mass onto the primary as a function 
of time for the 1.40+1.50 Mq binary with e = 0.15. The number of 
particles transferred during each episode is ~ 200, which represents 
an approximative lower limit to our ability to adequately resolve 
mass transfer due to the statistical nature of SPH. 



ergy is positive, the particle is unbound and is assigned 
to the ejecta. 

5.2. Estimates of mass transfer rates 

Using the total mass bound to the stellar components 
as a function of time, we can determine the mass transfer 
rates for the system modeled. We use a simple approach 
to determine the instantaneous mass transfer rates based 
on the difference of the total mass bound of each compo- 
nent between two successive timesteps, i.e. 



M 



Ml - Ml 
At 



(10) 



where M refers to the bound mass and the indices refer 
to component i and timesteps t and t — 1. 

Figure [10] shows the mass transfer rate of the primary 
as a function of time. Distinct episodes of mass transfer 
are observed to occur once per orbit and peak around 
a few 10 -6 Mq yr _1 . The number of particles trans- 
ferred during each episodes is 175 — 200 which, given the 
masses of the SPH particles, may limit our ability to re- 
solve lower mass transfer rates. Indeed, SPH requires a 
minimum number of neighbouring particles to calculate 
the density and in cases where only a handful of particles 
are transferred, the SPH treatment may not be adequate. 
Given the masses of the particles, low numbers of parti- 
cles therefore set lower limits to the mass transfer rates 
that our simulations can resolve. Using more particles of 
smaller masses would definitely allow for the resolution of 
lower mass transfer rates, as discussed in Sj6] Apart from 
the main episodes of mass transfer, we also observe sec- 
ondary peaks occurring before the main episodes of mass 
transfer. Our simulation shows that some of the mate- 
rial lost both by the primary and the secondary falls back 
onto both components and this is what is observed here. 



6. DISCUSSION 

Analytical approaches used to study the evolution of 
binaries usually rely on prescriptions or approximations 
when dealing with mass transfer rates. To relax some of 
these approximations, hydrodynamical modeling can be 
used, although it remains a hard task for many physi- 
cal and numerical reasons. To circumvent some of these 
difficulties, we introduced an alternate technique using 
boundary conditions and ghost particles to model only 
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the outermost parts of both the donor and the accre- 
tor stars. The location of the boundary is arbitrar- 
ily set. Since only the surface material is involved in 
mass transfer, our approach allows for better spatial and 
mass resolutions in the stream of matter. Moreover, our 
method allows for the modeling of both the accretor and 
the donor simultaneously while using less CPU time and 
maintaining realistic density profiles, taken from stellar 
models. 

Our code was shown to work particularly well for stars 
of equal mass and stars that are centrally condensed. In- 
deed, replacing the dense core of a massive star by a point 
mass is a good approximation. This is generally true for 
stars with masses > 0.8 M Q , although the evolutionary 
stage of the star may modify its density profile. Low- 
mass stars have shown to be more difficult to properly 
model with our approach (see Figure [7]) and we suggest 
that limiting our new boundary conditions to centrally 
condensed stars will in general yield better results. 

We also discussed the setup of proper initial conditions 
for modeling binary stars, which we have tested on both 
equal- and unequal-mass binaries. We demonstrated that 
our relaxation procedure is consistently implemented in 
our code and that it allows for the evolution over many 
orbits of equal-mass detached circular binaries and main- 
tain their orbital separation to within 1 — 2%. 

In light of the results from our first simulation of mass 
transfer, we establish that typical mass transfer rates 
that can be modeled with our new boundary conditions 
(i.e. > 10~ 6 Mp, yr" 1 ) a re consistent with the esti- 
mates of lChen &: Hani (|2008|) who investigated the forma- 
tion of blue stragglers through episodes of mass transfer 
onto main-sequence stars. We therefore suggest that our 
method consisting of modeling both stars simultaneously 
with appropriate boundary conditions can be applied to 
the problem of mass transfer in main-sequence binaries 
and help clarify the origin of blue stragglers. 

Using particles of lower mass allows for the modeling of 
lower mass transfer rates, although this requires the use 
of more particles and CPU time. Pushing the boundary 
furthe r out or using a poin t mass to model the secondary 
(as in IChurch et all I2009D can also allow for the use of 
more particles and a better mass transfer rate resolu- 
tion. Using a point mass as the secondary would however 
counter the benefit of our method to be able to model 
two interacting stars simultaneously. 

The physicality of our simulations depends of course 
on the physical ingredients we put in our code. As such, 
we do not include the effects of radiation pressure and 
energy transport mechanisms by radiation or convec- 
tion. These effects may have significance especially when 
studying the long-term evolution of mass-transferring bi- 
naries, where radiative cooling in the outer layers of the 
stars and envelope might be more important. Moreover, 
like any numerical technique, our method has some of its 
own limitations, and we discuss them now. 

6.1. Solid boundary 

By construction, our boundary is "semi-impermeable" , 
in that it does not allow particles to go through it. We set 
three smoothing lengths of ghosts and enforce that par- 
ticles be artificially repositioned if they happen to cross 
the boundary. However, it is possible that these con- 
ditions fail when dealing with large mass transfer rates 



and if some particles find themselves inside the bound- 
ary, around the central point mass, our code has to be 
stopped. This particle penetration limits our ability to 
adequately model episodes of extremely (and unrcalisti- 
cally) large mass transfer rates (> 10 _1 Mq yr^ 1 ; see 
Paper II). However, we do not think that our boundary 
should be so particle-tight since we do expect some mix- 
ing in the envelope of the secondaries. Although moving 
the boundary to a smaller radius could fix the issue of 
particle penetration, it would counter the use and ben- 
efits of our approach. Instead, we suggest using sink 
particles at the centre of the stars in order to account for 
deep mixing. Sink particles are like point masses but, 
in addition, their mass and momentum are allowed to 
increase as SPH particles get accreted. 

Also, as discussed in §3.2[ the relative position of our 
ghost particles is fixed in time, thus imitating a solid 
boundary. The boundary is not allowed to change its 
shape and/or provide a time- variable pressure gradient 
on the SPH particles. As a first approximation, thi s 
is a valid treatment (e.g. iDeupree fc Karakasl 120051) . 
However, when the gravitational potential changes sig- 
nificantly along the orbit, like on eccentric orbits, tidal 
forces may become non-negligible. But modeling the ef- 
fect of tidal forces on the boundary is costly, in terms 
of CPU time, as it involves calculating the gravitational 
force on the ghosts. This calculation is not done in our 
code as of now. We showed, however, that placing the 
boundary deep inside the star decreases the effect of 
tidal forces on the boundary and validates the use of 
our method. Finally, we note that the angular velocity 
of the ghost particles is maintained fixed during our sim- 
ulations. This is a valid assumption as synchronization 
occurs over timescales that are much longer (10 6 — 10 8 
years) than the duration of our simulations. 

6.2. Relaxation of unequal-mass binaries 

Our relaxation procedure for binary stars has proven to 
be especially efficient for equal-mass binaries. However, 
for unequal-mass binaries, we do not quite achieve the 
same level of accuracy for the evolution of orbital sep- 
aration. We think the reason for this difference comes 
from the fact that the equal-mass systems we model are 
perfectly symmetric, i.e. the two stars are exact replicas 
of each other, whereas in the case of unequal-mass bina- 
ries, symmetry is broken. For equal-mass binaries, the 
gravitational acceleration calculations are exactly equal 
and opposite and the two stars are evolved identically. 
But given the adaptive nature of our code, stars (and 
particles) of different masses may be evolved on differ- 
ent timesteps and care should be taken if the two stars 
(and particles) are to be evolved consistently. We per- 
formed test runs during which we forced our code to use 
a smaller common timestep, but this approach does not 
improve the results. On the other hand, using a more di- 
rect summation approach for the calculation of the grav- 
itational force (by using a smaller binary tree opening 
angle) is found to improve the results. Indeed, results 
from test runs using such an approach show much im- 
provements in evolving stars on circular orbits. However, 
doing so also makes our simulations significantly longer 
to run. Therefore, we think the observed behaviour of 
our unequal-mass binaries is the result of our calcula- 
tion of gravitational acceleration through a binary tree, 
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although more work remains to be done. 

6.3. Future Work 

The method presented in this paper has been shown 
to be well suited for modeling the hydrodynamics of in- 
teracting binary stars. Here, we used it to model Roche 
lobe overflow, during which only the outermost layers 
of the stars are actively involved. We propose that the 
approach presented in this work can also be applied to 
many different situations, such as wind accretion and in- 
terstellar medium accretion, and that it can help better 
understand how stars react to mass loss and mass accre- 
tion in general. In particular, we have emphasized that 
the Roche lobe formalism is not applicable in the case 
of eccentric and asynchronous binaries. We think our 



alternate method can be used to better understand and 
characterize the onset of mass transfer in such systems. 
This is the subject of a subsequent paper (jLaioie fc Sills) 
l2010f ) in which we model eccentric systems of different 
masses, semi-major axes, and eccentricities. 
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